1 Présentation
2 Import des données
3 Préparation des données
4 Étude auto-corrélation et colinéarité
- 4.1 Réalisation d’une matrice de corrélation
- 4.2 Calcul de la Variance Inflaction Factor (VIF)
5 Comparaison des valeurs des communes fusionnantes vis-à-vis des valeurs des communes inchangées
6 Régression logistique classique
7 [À faire] Étude des résidus (y compris avec indice de Moran pour une auto-corrélation spatiale des résidus ?)
8 [À vérifier, tester] Forward stepwise ou backward stepwise
9 Bibliographie

1 Présentation

L’objectif de ce script est de tester différents modèles de régressions logistiques permettant de modéliser le phénomène des communes nouvelles.

Ce travail s’inscrit dans le cadre de la thèse de G. Bideau sur les communes nouvelles, sous la direction de R. Le Goix.

Plan envisagé au départ :

Import de différentes données (par blocs complets)
Prévoir des tests sur des sous-ensembles (Maine-et-Loire, Savoies, Ouest…)Il est important voire intéressant de cartographier ou de représenter graphiquement les résidus. Cela permet potentiellement de valider ou invalider le modèle mais également de réfléchir aux espaces où ce modèle fonctionne le mieux ou le moins bien.
Calculs colinéarité, corrélations
Méthode de régression classique (copier-coller analyse budgets)
Forward stepwise ou backward stepwise
Cartographie des résidus (y compris avec indice de Moran pour une auto-corrélation spatiale des résidus ?)

2 Import des données

2.1 Géométries

On n’importe au départ que les données à la géométrie 2011, mentionnant les futures communes nouvelles de rattachement.

# Import des géométries 2011
geom2011 <- st_read("Archives/data_prep 2011-2020(01)/data/geom.gpkg", layer = "geom2011", quiet = TRUE) 
geomfus2011 <- st_read("Archives/data_prep 2011-2020(01)/data/geom.gpkg", layer = "geomfus2011", quiet = TRUE) 
dep <- st_read("Archives/data_prep 2011-2020(01)/data/geom.gpkg", layer = "dep", quiet = TRUE)
geom_new <- st_read("Archives/data_prep 2011-2020(01)/data/geom.gpkg", layer = "geom_new", quiet = TRUE) 
geomCN_new <- st_read("Archives/data_prep 2011-2020(01)/data/geom.gpkg", layer = "geomCN_new", quiet = TRUE)

2.2 Base DAC

Base DAC [bideau2022b]

# Import des données de la base DAC
load("data/refdata.Rdata")

variables_RT <- colnames(df2011[which(str_detect(colnames(df2011), "_RT") == TRUE)])

2.3 Budgets

# Import des données budgétaires pour l'année 2011
load("data/refdata_budgets_bloccommunal_tot.Rdata")
rm(df_budgets_new, budgets_com_2011, budgets_com_2020)

# Fusion des données budgétaires avec les autres données
# On liste les variables renseignées dans les différents df
colbudgets2011 <- colnames(df_budgets_2011)

coldf2011 <- colnames(df2011)

# On liste les variables présentes seulement dans budgets (celles qu'il va falloir y prélever)
colgarder2011 <- setdiff(colbudgets2011, coldf2011)

# On fusionne les données de data_prep et des budgets, en ne gardant pour ces dernières que les données précédemment identifiées
df2011 <- merge(df2011,
                     df_budgets_2011[, c("CODGEO", colgarder2011)],
                     by = "CODGEO", all.x = TRUE, all.y = TRUE)

# Toutes les variables en pourcentages qui sont des variables budgétaires
variables_budg_prct <- colnames(df2011[which(str_detect(colnames(df2011), "_prct") == TRUE)])
variables_budg_prct <- variables_budg_prct[variables_budg_prct!= "geometry"] # On supprime la variable "geometry" qui est intégrée sans que cela soit voulu


# Suppression des données inutiles (déjà importées ou seulement de transition)
rm(df_budgets_2011, colbudgets2011, colgarder2011)

2.4 Données électorales 2012

PR2012_T1 <- read.table("data/elections/PR2012_T1.csv", sep="\t", colClasses = c(rep("character", 2), rep("numeric", 40)), head = TRUE, stringsAsFactors = TRUE, dec =",")
colnames(PR2012_T1)[3:42] <- paste0("PR2012_T1_", colnames(PR2012_T1)[3:42])

# PR2012_T2 <- read.table("data/elections/PR2012_T2.csv", sep="\t", colClasses = c(rep("character", 2), rep("numeric", 16)), head = TRUE, stringsAsFactors = TRUE, dec =",")

variables_elect <- c("PR2012_T1_Abst_prct_insc", "PR2012_T1_prct_insc_LE.PEN", "PR2012_T1_prct_insc_SARKOZY", "PR2012_T1_prct_insc_MÉLENCHON", "PR2012_T1_prct_insc_BAYROU", "PR2012_T1_prct_insc_HOLLANDE")

colPR2012_T1 <- colnames(PR2012_T1)
# colPR2012_T2 <- colnames(PR2012_T2)

# On liste les variables présentes seulement dans budgets (celles qu'il va falloir y prélever)
colgarder2011_PR2012_T1 <- setdiff(colPR2012_T1, coldf2011)
# colgarder2011_PR2012_T2 <- setdiff(colPR2012_T2, coldf2011)

# On fusionne les données de data_prep et des budgets, en ne gardant pour ces dernières que les données précédemment identifiées
df2011 <- merge(df2011,
                     PR2012_T1[, c("CODGEO", colgarder2011_PR2012_T1)],
                     by = "CODGEO", all.x = TRUE, all.y = TRUE)
# Suppression des données inutiles (déjà importées ou seulement de transition)
# df2011 <- merge(df2011,
#                      PR2012_T2[, c("CODGEO", colgarder2011_PR2012_T2)],
#                      by = "CODGEO", all.x = TRUE, all.y = TRUE)
rm(PR2012_T1, PR2012_T2, colPR2012_T1, colPR2012_T2, coldf2011, colgarder2011_PR2012_T1, colgarder2011_PR2012_T2)

## Warning in rm(PR2012_T1, PR2012_T2, colPR2012_T1, colPR2012_T2, coldf2011, :
## objet 'PR2012_T2' introuvable

## Warning in rm(PR2012_T1, PR2012_T2, colPR2012_T1, colPR2012_T2, coldf2011, :
## objet 'colPR2012_T2' introuvable

## Warning in rm(PR2012_T1, PR2012_T2, colPR2012_T1, colPR2012_T2, coldf2011, :
## objet 'colgarder2011_PR2012_T2' introuvable

3 Préparation des données

3.1 Sélection études de cas

Section qui sert au départ à tester le code sur de petits ensembles.

À l’avenir, pourra servir à définir des boucles

EdCs <- c("EdC_France", "EdC_Savoies", "EdC_Ouest", "EdC_Normandie", "EdC_49", "EdC_RALP_partiel")

df2011$EdC_France <- "OUI" # France entière

df2011$EdC_Savoies[df2011$CODE_DEPT %in% c("73", "74")] <- "OUI" # Départements de Savoie et Haute-Savoie
df2011[, c("EdC_Savoies")] [is.na(df2011[, c("EdC_Savoies")])] <- "NON"

df2011$EdC_Ouest[df2011$REG %in% c("23", "25", "53", "52")] <- "OUI" # Haute-Normandie, Basse Normandie, Bretagne, Pays-de-la-Loire
df2011[, c("EdC_Ouest")] [is.na(df2011[, c("EdC_Ouest")])] <- "NON"

df2011$EdC_Normandie[df2011$REG %in% c("23", "25")] <- "OUI" # Haute et Basse Normandie
df2011[, c("EdC_Normandie")] [is.na(df2011[, c("EdC_Normandie")])] <- "NON"

df2011$EdC_49[df2011$CODE_DEPT %in% c("49")] <- "OUI" # Département du Maine-et-Loire
df2011[, c("EdC_49")] [is.na(df2011[, c("EdC_49")])] <- "NON"

df2011$EdC_RALP_partiel[df2011$CODE_DEPT %in% c("69", "73", "74", "38", "01")] <- "OUI" # Départements de Savoie, Haute-Savoie, Isère, Rhône et Ain
df2011[, c("EdC_RALP_partiel")] [is.na(df2011[, c("EdC_RALP_partiel")])] <- "NON"

EdC <- EdCs[1]

df_reg <- subset(df2011, df2011[EdC] == "OUI")
df_reg_Cfus <- subset(df_reg, df_reg$COM_NOUV == "OUI")
# On peut choisir de se limiter aux communes françaises ayant une population ne dépassant pas celle de la plus peuplée des communes nouvelles
df_reg <- subset(df_reg, df_reg$P09_POP <= max(df_reg_Cfus$P09_POP))

# On ajoute les géométries
geomfus2011 <- merge(geomfus2011, df_reg[, c("CODGEO", "LIBGEO")], by = "CODGEO", all.x = FALSE, all.y = TRUE)
df_reg <- merge(geom2011, df_reg, by = "CODGEO", all.y = TRUE)



plot(df_reg$geom, main = paste("Espace étudié :", EdC))

3.2 Adaptation variables pourcentages

NB : Il y a un problème si on utilise les variables budgétaires en valeurs absolues ou les variables en pourcentages notées 100% = 100. Cela fonctionne correctement si on a les pourcentages notées en 100% = 1.

Deux possibilités : diviser les pourcentages manuellement ou utiliser le fichier ‘refdata_budgets_prct_en_ratio.Rdata’. On choisit la première solution.

Pas vraiment nécessaire car passage des variables quantitatives en variables qualitatives (quantiles).

# Pour multiplier/diviser
# Toutes les variables en pourcentages qui sont des variables budgétaires
variables_budg_prct <- colnames(df_reg[which(str_detect(colnames(df_reg), "_prct") == TRUE)])
variables_budg_prct <- variables_budg_prct[variables_budg_prct!= "geometry"] # On supprime la variable "geometry" qui est intégrée sans que cela soit voulu
df_reg [variables_budg_prct] <- df_reg[variables_budg_prct]/100
# test <- df_reg[variables_budg_prct]

variables_RT <- colnames(df_reg[which(str_detect(colnames(df_reg), "_RT") == TRUE)])
variables_RT <- variables_RT[variables_RT!= "geometry"] # On supprime la variable "geometry" qui est intégrée sans que cela soit voulu
df_reg [variables_RT] <- df_reg[variables_RT]/100
# test <- df_reg[variables_RT]

# Sans que je ne comprenne bien pourquoi, les géométries sont également touchées, si variable non retirée avant, on rétablit donc ces dernières
# df_reg$geometry <- df_reg$geometry * 10000

3.3 Transformations variables quantitatives en variables qualitatives

Le principe d’une régression logistique est de définir la probabilité d’un évènement (dans mon cas, fusionner ou non) en fonction du changement d’une variable. Dans le cas d’une variable quantitative, la question est quelle est le pourcentage de chance que l’évènement survienne davantage lorsqu’on augmente la variable quantitative de 1.

Dans le cas d’une variable qualitative, il s’agit du pourcentage de chance que l’évènement survienne davantage lorsqu’on change de catégorie. Il faut donc faire attention à mettre la variable en facteur, avec la première valeur qui doit être celle de référence (c’est à partir de celle-là que les comparaisons vont être faites, c’est donc important car c’est cela qui détermine le discours).

# variable <- "C09_ACTOCC_OUT_RT"
variables_quanti_quali <- c(variables_budg_prct, variables_RT, variables_elect)
# On retire certaines variables qui ne peuvent être divisées par tertiles car trop de valeurs nulles.
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DSU_tot_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DNP_princip_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DNP_sortie_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DNP_major_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DNP_tot_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DSR_bourg_centre_global_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="DSR_cible_global_prct"]
variables_quanti_quali <- variables_quanti_quali[variables_quanti_quali!="C09_EMPLT_INDUS_RT"]

df_sans_geom <- df_reg # Création d'un jeu de données pour le travail
st_geometry(df_sans_geom) <- NULL # Suppression des géométries
# summary(df_sans_geom[, variable])

for (variable in variables_quanti_quali) {

Y <- df_sans_geom[, variable] # On se focalise sur une variable en particulier

# Calcul moyenne groupes
moy <- as.data.frame(tapply(df_sans_geom[, variable], df_sans_geom$COM_NOUV, mean, na.rm = TRUE))
moy$COM_NOUV <- row.names(moy)
colnames(moy) <- c("moyenne", "COM_NOUV") 

# Calcul moyenne groupes
med <- as.data.frame(tapply(df_sans_geom[, variable], df_sans_geom$COM_NOUV, median, na.rm = TRUE))
med$COM_NOUV <- row.names(moy)
colnames(med) <- c("median", "COM_NOUV") 

# On veut observer graphiquement la distribution statistique
histo <- ggplot(df_sans_geom)+
  geom_histogram(aes(x=df_sans_geom[, variable], color = COM_NOUV), bins=50)+
  labs(title = paste("Variable étudiée :\n", variable), subtitle = "Pointillés : moyennes.\nLigne pleine : médianes.")+
  # geom_vline(aes(xintercept=median(df_sans_geom[, variable], na.rm = T)),      color="black", linetype="dashed" )+
  geom_vline(data=moy, aes(xintercept=moyenne, color=COM_NOUV),             linetype="dashed") +
  geom_vline(data=med, aes(xintercept=median, color=COM_NOUV)) +
  ylab("Count")+theme_light()
print(histo)

# On découpe en fonction du quantile souhaité puis on modifie l'ordre de niveaux de facteur pour que la variable centrale soit la première

coupures <- quantile(Y, probs=seq(0, 1, 0.2), na.rm = TRUE)

if (any(duplicated(coupures)) == FALSE) { # Si c'est possible, on utilise les quintiles
  Y <- cut(Y,breaks=c(quantile(Y, probs=seq(0, 1, 0.2), na.rm = TRUE))) # Pour découpage en quintiles
levels(Y)<-c("Q1","Q2", "Q3", "Q4", "Q5")
Y <- relevel (Y, "Q3","Q1", "Q2", "Q4", "Q5")
}else{# Sinon, on utilise les tertiles
Y <- cut(Y,breaks=c(quantile(Y, probs=seq(0, 1, 1/3), na.rm = TRUE))) # Pour découpage en tertiles
levels(Y)<-c("Tertile_inf","Tertile_med", "Tertile_sup")
Y <- relevel (Y, "Tertile_med", "Tertile_inf", "Tertile_sup")
}


# Y<-cut(Y,breaks=c(quantile(Y))) # Pour découpage en quartiles
# levels(Y)<-c("Q1","Q2", "Q3", "Q4")




summary(Y)


X <- df_sans_geom$COM_NOUV
summary(X)

tabcont<-table(X,Y)
tabcont # En valeur absolue
# round(100*prop.table(tabcont,margin=1),1) # Pourcentages, le total se fait par lignes
# round(100*prop.table(tabcont,margin=),1) # Pourcentages, le total se fait sur l'ensemble de la population
# round(100*prop.table(tabcont,margin=2),1) # Pourcentages, le total se fait par colonnes

# On verse les nouvelles données au data frame
df_reg[, paste0(variable, "_quali")] <- Y

}

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 11 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 52 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 52 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 52 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 101 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 126 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 958 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_bin()`).

rm(tabcont, X, Y, df_sans_geom, variable, histo)

3.4 Modifications de quelques autres variables

# On réordonne le facteur CATAEU2010
df_reg$CATAEU2010 <- as.factor(df_reg$CATAEU2010)
freq(df_reg$CATAEU2010, sort = "dec")

##         n    % val%
## 112 12180 33.7 33.7
## 400  7233 20.0 20.0
## 300  6975 19.3 19.3
## 120  3968 11.0 11.0
## 111  3100  8.6  8.6
## 221   854  2.4  2.4
## 212   801  2.2  2.2
## 222   558  1.5  1.5
## 211   439  1.2  1.2

row.names(freq(df_reg$CATAEU2010, sort = "dec"))

## [1] "112" "400" "300" "120" "111" "221" "212" "222" "211"

# Vérifier que la modalité qui se trouve en premier est bien la plus fréquente.
# Sinon, possibilité d'utiliser le code suivant :
# df_reg$CATAEU2010 <- relevel (df_reg$CATAEU2010, "112", "111", "120", "211", "212", "221", "222", "300", "400")
# df_reg$CATAEU2010 <- relevel (df_reg$CATAEU2010, "112", "400", "300", "120", "111", "221", "212", "222", "211")
ordre <- row.names(freq(df_reg$CATAEU2010, sort = "dec"))
df_reg$CATAEU2010 <- factor (df_reg$CATAEU2010, ordre)


df_reg$CODE_DEPT <- as.factor(df_reg$CODE_DEPT)
freq(df_reg$CODE_DEPT, sort = "dec")

##      n   % val%
## 62 894 2.5  2.5
## 02 815 2.3  2.3
## 80 781 2.2  2.2
## 76 743 2.1  2.1
## 57 729 2.0  2.0
## 14 705 2.0  2.0
## 21 705 2.0  2.0
## 60 692 1.9  1.9
## 27 675 1.9  1.9
## 59 645 1.8  1.8
## 51 619 1.7  1.7
## 50 601 1.7  1.7
## 25 593 1.6  1.6
## 54 593 1.6  1.6
## 31 588 1.6  1.6
## 71 573 1.6  1.6
## 24 557 1.5  1.5
## 64 546 1.5  1.5
## 70 545 1.5  1.5
## 39 544 1.5  1.5
## 33 539 1.5  1.5
## 38 532 1.5  1.5
## 67 526 1.5  1.5
## 88 515 1.4  1.4
## 77 513 1.4  1.4
## 61 505 1.4  1.4
## 55 500 1.4  1.4
## 65 474 1.3  1.3
## 17 471 1.3  1.3
## 63 469 1.3  1.3
## 08 463 1.3  1.3
## 32 463 1.3  1.3
## 89 455 1.3  1.3
##  [ reached 'max' / getOption("max.print") -- omitted 60 rows ]

ordre <- row.names(freq(df_reg$CODE_DEPT, sort = "dec"))
df_reg$CODE_DEPT <- factor (df_reg$CODE_DEPT, ordre)

df_reg$REG <- as.factor(df_reg$REG)
freq(df_reg$REG, sort = "dec")

##       n   % val%
## 73 3018 8.4  8.4
## 82 2872 8.0  8.0
## 41 2337 6.5  6.5
## 72 2292 6.3  6.3
## 22 2288 6.3  6.3
## 26 2045 5.7  5.7
## 21 1947 5.4  5.4
## 24 1839 5.1  5.1
## 25 1811 5.0  5.0
## 43 1784 4.9  4.9
## 91 1541 4.3  4.3
## 31 1539 4.3  4.3
## 52 1497 4.1  4.1
## 54 1459 4.0  4.0
## 23 1418 3.9  3.9
## 83 1309 3.6  3.6
## 53 1265 3.5  3.5
## 11 1247 3.5  3.5
## 93  953 2.6  2.6
## 42  901 2.5  2.5
## 74  746 2.1  2.1

row.names(freq(df_reg$REG, sort = "dec"))

##  [1] "73" "82" "41" "72" "22" "26" "21" "24" "25" "43" "91" "31" "52" "54" "23"
## [16] "83" "53" "11" "93" "42" "74"

ordre <- row.names(freq(df_reg$REG, sort = "dec"))
df_reg$REG <- factor (df_reg$REG, ordre)


df_reg$P09_POP <- as.numeric(df_reg$P09_POP)
df_reg$P09_POP_quali <- cut(df_reg$P09_POP, breaks=c(quantile(df_reg$P09_POP, probs=seq(0, 1, 1/5), na.rm = TRUE)))
levels(df_reg$P09_POP_quali)<-c("Q1","Q2", "Q3", "Q4", "Q5")
df_reg$P09_POP_quali <- relevel (df_reg$P09_POP_quali, "Q3","Q1", "Q2", "Q4", "Q5")

4 Étude auto-corrélation et colinéarité

Pour éviter la colinéarité entre les variables explicatives, deux méthodes [feuillet2019] :

faire une matrice de corrélation et supprimer les variables dont les corrélations sont supérieures à 0,7 en valeur absolue [tenenhaus2007] ;
calculer la Variance Inflation Factor et supprimer les variables avec une VIF supérieure à 3 (voire supérieur à 2, même référence biblio).

4.1 Réalisation d’une matrice de corrélation

Cf. https://cran.r-project.org/web/packages/corrplot/vignettes/corrplot-intro.html

# Pour sélectionner certaines variables spécifiquement
# variables_correlations <- c("fctva_prct", "dgf_prct", "DSU_tot_prct", "DNP_tot_prct", "DSR_pereq_global_prct")
# variables_correlations <- c("fctva", "dgf", "DSU_tot", "DNP_tot", "DSR_pereq_global")
# variables_correlations <- vartests

# Toutes les variables en pourcentages NB : ajout du symbole $ pour indiquer que c'est la fin de la chaîne de caractère
variables_correlations <- c(colnames(df_reg[which(str_detect(colnames(df_reg), "_prct$") == TRUE)]),
                            colnames(df_reg[which(str_detect(colnames(df_reg), "_RT$") == TRUE)]))
variables_correlations <- c(colnames(df_reg[which(str_detect(colnames(df_reg), "_prct$") == TRUE)]),
                            colnames(df_reg[which(str_detect(colnames(df_reg), "_RT$") == TRUE)]), variables_elect)

variables_correlations <- variables_correlations[variables_correlations!= "geometry"] # On supprime la variable "geometry" si elle est intégrée sans que cela soit voulu

# variables_correlations <- c("P09_RETR1564_RT", "C09_ACT1564_Agr_RT", "P09_POP3044Y_RT", "C09_ACTOCC_OUT_RT", "P09_CHOM1564_RT", "dgf_prct", "perso_prct", "charge_prct")

# Définition du seuil d'auto-corrélation
seuil <- 0.7

# Juste pour les communes fusionnées
tmp <- subset(df2011, df2011$COM_NOUV == "OUI")
tmp <- na.omit(tmp [, variables_correlations])

tmp2 <- cor(tmp, method = "pearson")
res1 <- cor.mtest(tmp, conf.level = .99)

col <- c("#a50026","#d73027","#f46d43","#fdae61","#fee090","#e0f3f8","#abd9e9","#74add1","#4575b4","#313695")

corrplot::corrplot(tmp2, p.mat = res1$p, sig.level = .2, order = "hclust", tl.cex = 0.8, tl.col = "black") # Obligé d'appeler le package directement car sinon conflit avec une autre fonction du package `arm`.

corrplot::corrplot(tmp2, method="color", order = "hclust", # cl.lim = c(-1,1), # cl.lim pose problème donc supprimé
         col = col,
         tl.col = rgb(0,0,0),
         number.cex = 0.4,
         addCoef.col = rgb(0,0,0, alpha =0.6),
         addgrid.col = "white",
         title = "Analyse des corrélations",
         tl.cex=0.8,
         cl.cex = 0.6,
         mar=c(0,0,4,0))

tmp2[which(tmp2 == 1)] <- 0 # On remplace par 0 l'autocorrélations des variables avec elles-mêmes pour faciliter la lecture ensuite
# which(abs(tmp2) > 0.7, arr.ind=TRUE) # On identifie les valeurs dont la valeur absolue serait supérieure à 0,7

for(i in 1:nrow(tmp2)) {
   for(j in 1:ncol(tmp2))
   {if(tmp2[i,j]>seuil)
   {print(paste("Variables auto-corrélées à plus de", seuil, ":", rownames(tmp2)[i], "et", colnames(tmp2)[j], ":", tmp2[i,j]))
     }}}

## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et charge_prct : 0.732689123376272"
## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et res1_prct : 0.734728134238698"
## [1] "Variables auto-corrélées à plus de 0.7 : impo1_prct et bfb_prct : 0.934797223269612"
## [1] "Variables auto-corrélées à plus de 0.7 : impo1_prct et pfb_prct : 0.944226808552536"
## [1] "Variables auto-corrélées à plus de 0.7 : impo1_prct et PotFin4taxes_prct : 0.723731388000992"
## [1] "Variables auto-corrélées à plus de 0.7 : impo1_prct et P11_POT_FIN_RT : 0.823753548373962"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DotationNmoins1PerimN_prct : 0.958902817093769"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et Dotation_forfaitaireN_prct : 0.96347828612023"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DSR_pereq_global_prct : 0.826044221470363"
## [1] "Variables auto-corrélées à plus de 0.7 : charge_prct et prod_prct : 0.732689123376272"
## [1] "Variables auto-corrélées à plus de 0.7 : res1_prct et prod_prct : 0.734728134238698"
## [1] "Variables auto-corrélées à plus de 0.7 : recinv_prct et res2_prct : 0.741258939982701"
## [1] "Variables auto-corrélées à plus de 0.7 : depinv_prct et equip_prct : 0.887426789236825"
## [1] "Variables auto-corrélées à plus de 0.7 : equip_prct et depinv_prct : 0.887426789236825"
## [1] "Variables auto-corrélées à plus de 0.7 : remb_prct et annu_prct : 0.97130987424642"
## [1] "Variables auto-corrélées à plus de 0.7 : res2_prct et recinv_prct : 0.741258939982701"
## [1] "Variables auto-corrélées à plus de 0.7 : annu_prct et remb_prct : 0.97130987424642"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et impo1_prct : 0.934797223269612"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et pfb_prct : 0.987194884577028"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et PotFin4taxes_prct : 0.722401553651033"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et P11_POT_FIN_RT : 0.87843217078898"
## [1] "Variables auto-corrélées à plus de 0.7 : bfnb_prct et pfnb_prct : 0.736196383226621"
## [1] "Variables auto-corrélées à plus de 0.7 : pfb_prct et impo1_prct : 0.944226808552536"
## [1] "Variables auto-corrélées à plus de 0.7 : pfb_prct et bfb_prct : 0.987194884577028"
## [1] "Variables auto-corrélées à plus de 0.7 : pfb_prct et P11_POT_FIN_RT : 0.884623222219747"
## [1] "Variables auto-corrélées à plus de 0.7 : pfnb_prct et bfnb_prct : 0.736196383226621"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et dgf_prct : 0.958902817093769"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et Dotation_forfaitaireN_prct : 0.996833840765774"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et DSR_pereq_global_prct : 0.804191406661557"
## [1] "Variables auto-corrélées à plus de 0.7 : PotFin4taxes_prct et impo1_prct : 0.723731388000992"
## [1] "Variables auto-corrélées à plus de 0.7 : PotFin4taxes_prct et bfb_prct : 0.722401553651033"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et dgf_prct : 0.96347828612023"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DotationNmoins1PerimN_prct : 0.996833840765774"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DSR_pereq_global_prct : 0.812335469877741"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_major_prct : 0.780156401039296"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_tot_prct : 0.97878473516058"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_princip_prct : 0.780156401039296"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_tot_prct : 0.890785442109139"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_princip_prct : 0.97878473516058"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_major_prct : 0.890785442109139"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et dgf_prct : 0.826044221470363"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et DotationNmoins1PerimN_prct : 0.804191406661557"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et Dotation_forfaitaireN_prct : 0.812335469877741"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_RETR1564_RT et P09_POP6074Y_RT : 0.742995991320959"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_ACT1564_Agr_RT et C09_EMPLT_AGRI_RT : 0.722572663768533"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_EMPLT_AGRI_RT et C09_ACT1564_Agr_RT : 0.722572663768533"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_POP6074Y_RT et P09_RETR1564_RT : 0.742995991320959"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_POT_FIN_RT et impo1_prct : 0.823753548373962"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_POT_FIN_RT et bfb_prct : 0.87843217078898"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_POT_FIN_RT et pfb_prct : 0.884623222219747"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_IMP_NET_RT : 0.861018323277089"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_FoyFisc_Imp_RT : 0.738703335772406"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_IMP_NET_RT et P11_Rev_Fisc_RT : 0.861018323277089"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_FoyFisc_Imp_RT et P11_Rev_Fisc_RT : 0.738703335772406"

# Communes n'ayant pas fusionné
tmp <- subset(df2011, df2011$COM_NOUV == "NON")
tmp <- na.omit(tmp [, variables_correlations])

tmp2 <- cor(tmp, method = "pearson")
res1 <- cor.mtest(tmp, method = "pearson")

corrplot::corrplot(tmp2, p.mat = res1$p, sig.level = .2, order = "hclust",
         tl.cex = 0.8, tl.col = "black")

col <- c("#a50026","#d73027","#f46d43","#fdae61","#fee090","#e0f3f8","#abd9e9","#74add1","#4575b4","#313695")

corrplot::corrplot(tmp2, method="color", order = "hclust", # cl.lim = c(-1,1), # cl.lim pose problème donc supprimé 
         col = col,
         tl.col = rgb(0,0,0),
         number.cex = 0.4,
         addCoef.col = rgb(0,0,0, alpha =0.6),
         addgrid.col = "white",
         title = "Analyse des corrélations",
         tl.cex=0.8,
         cl.cex = 0.6,
         mar=c(0,0,4,0))

tmp2[which(tmp2 == 1)] <- 0 # On remplace par 0 l'autocorrélations des variables avec elles-mêmes pour faciliter la lecture ensuite
# which(abs(tmp2) > 0.7, arr.ind=TRUE) # On identifie les valeurs dont la valeur absolue serait supérieure à 0,7

for(i in 1:nrow(tmp2)) {
   for(j in 1:ncol(tmp2))
   {if(tmp2[i,j]>seuil)
   {print(paste("Variables auto-corrélées à plus de", seuil, ":", rownames(tmp2)[i], "et", colnames(tmp2)[j], ":", tmp2[i,j]))
     }}}

## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et charge_prct : 0.724234792513559"
## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et res1_prct : 0.732830571012828"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DotationNmoins1PerimN_prct : 0.945250511144103"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et Dotation_forfaitaireN_prct : 0.953487917133836"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DSR_pereq_global_prct : 0.79232599936631"
## [1] "Variables auto-corrélées à plus de 0.7 : charge_prct et prod_prct : 0.724234792513559"
## [1] "Variables auto-corrélées à plus de 0.7 : res1_prct et prod_prct : 0.732830571012828"
## [1] "Variables auto-corrélées à plus de 0.7 : recinv_prct et res2_prct : 0.723036360221147"
## [1] "Variables auto-corrélées à plus de 0.7 : depinv_prct et equip_prct : 0.902821742043709"
## [1] "Variables auto-corrélées à plus de 0.7 : equip_prct et depinv_prct : 0.902821742043709"
## [1] "Variables auto-corrélées à plus de 0.7 : remb_prct et annu_prct : 0.971637705415465"
## [1] "Variables auto-corrélées à plus de 0.7 : res2_prct et recinv_prct : 0.723036360221147"
## [1] "Variables auto-corrélées à plus de 0.7 : annu_prct et remb_prct : 0.971637705415465"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et PotFin4taxes_prct : 0.782069994452378"
## [1] "Variables auto-corrélées à plus de 0.7 : bfnb_prct et pfnb_prct : 0.733021436204166"
## [1] "Variables auto-corrélées à plus de 0.7 : pfnb_prct et bfnb_prct : 0.733021436204166"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et dgf_prct : 0.945250511144103"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et Dotation_forfaitaireN_prct : 0.991260468887141"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et DSR_pereq_global_prct : 0.771690289546452"
## [1] "Variables auto-corrélées à plus de 0.7 : PotFin4taxes_prct et bfb_prct : 0.782069994452378"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et dgf_prct : 0.953487917133836"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DotationNmoins1PerimN_prct : 0.991260468887141"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DSR_pereq_global_prct : 0.777669035129182"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_major_prct : 0.79425679995476"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_tot_prct : 0.979844133560585"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_princip_prct : 0.79425679995476"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_tot_prct : 0.897902584225198"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_princip_prct : 0.979844133560585"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_major_prct : 0.897902584225198"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et dgf_prct : 0.79232599936631"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et DotationNmoins1PerimN_prct : 0.771690289546452"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et Dotation_forfaitaireN_prct : 0.777669035129182"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_RETR1564_RT et P09_POP6074Y_RT : 0.740842384941785"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_ACT1564_Agr_RT et C09_EMPLT_AGRI_RT : 0.713423404098747"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_EMPLT_AGRI_RT et C09_ACT1564_Agr_RT : 0.713423404098747"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_POP6074Y_RT et P09_RETR1564_RT : 0.740842384941785"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_IMP_NET_RT : 0.880125493065463"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_FoyFisc_Imp_RT : 0.734113201668045"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_IMP_NET_RT et P11_Rev_Fisc_RT : 0.880125493065463"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_FoyFisc_Imp_RT et P11_Rev_Fisc_RT : 0.734113201668045"

# Ensemble des communes françaises
tmp <- df2011
tmp <- na.omit(tmp [, variables_correlations])

tmp2 <- cor(tmp, method = "pearson")
res1 <- cor.mtest(tmp, method = "pearson")

corrplot::corrplot(tmp2, p.mat = res1$p, sig.level = .2, order = "hclust",
         tl.cex = 0.8, tl.col = "black")

col <- c("#a50026","#d73027","#f46d43","#fdae61","#fee090","#e0f3f8","#abd9e9","#74add1","#4575b4","#313695")

corrplot::corrplot(tmp2, method="color", order = "hclust", # cl.lim = c(-1,1), # cl.lim pose problème donc supprimé 
         col = col,
         tl.col = rgb(0,0,0),
         number.cex = 0.4,
         addCoef.col = rgb(0,0,0, alpha =0.6),
         addgrid.col = "white",
         title = "Analyse des corrélations",
         tl.cex=0.8,
         cl.cex = 0.6,
         mar=c(0,0,4,0))

tmp2[which(tmp2 == 1)] <- 0 # On remplace par 0 l'autocorrélations des variables avec elles-mêmes pour faciliter la lecture ensuite
# which(abs(tmp2) > 0.7, arr.ind=TRUE) # On identifie les valeurs dont la valeur absolue serait supérieure à 0,7

for(i in 1:nrow(tmp2)) {
   for(j in 1:ncol(tmp2))
   {if(tmp2[i,j]>seuil)
   {print(paste("Variables auto-corrélées à plus de", seuil, ":", rownames(tmp2)[i], "et", colnames(tmp2)[j], ":", tmp2[i,j]))
     }}}

## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et charge_prct : 0.724818656916716"
## [1] "Variables auto-corrélées à plus de 0.7 : prod_prct et res1_prct : 0.732938988586129"
## [1] "Variables auto-corrélées à plus de 0.7 : impo1_prct et pfb_prct : 0.719361456090989"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DotationNmoins1PerimN_prct : 0.946561903962206"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et Dotation_forfaitaireN_prct : 0.954457043664577"
## [1] "Variables auto-corrélées à plus de 0.7 : dgf_prct et DSR_pereq_global_prct : 0.795902189919594"
## [1] "Variables auto-corrélées à plus de 0.7 : charge_prct et prod_prct : 0.724818656916716"
## [1] "Variables auto-corrélées à plus de 0.7 : res1_prct et prod_prct : 0.732938988586129"
## [1] "Variables auto-corrélées à plus de 0.7 : recinv_prct et res2_prct : 0.724291657635224"
## [1] "Variables auto-corrélées à plus de 0.7 : depinv_prct et equip_prct : 0.901768326955515"
## [1] "Variables auto-corrélées à plus de 0.7 : equip_prct et depinv_prct : 0.901768326955515"
## [1] "Variables auto-corrélées à plus de 0.7 : remb_prct et annu_prct : 0.971615494510789"
## [1] "Variables auto-corrélées à plus de 0.7 : res2_prct et recinv_prct : 0.724291657635224"
## [1] "Variables auto-corrélées à plus de 0.7 : annu_prct et remb_prct : 0.971615494510789"
## [1] "Variables auto-corrélées à plus de 0.7 : bfb_prct et pfb_prct : 0.878040662678628"
## [1] "Variables auto-corrélées à plus de 0.7 : bfnb_prct et pfnb_prct : 0.734066280198133"
## [1] "Variables auto-corrélées à plus de 0.7 : pfb_prct et impo1_prct : 0.719361456090989"
## [1] "Variables auto-corrélées à plus de 0.7 : pfb_prct et bfb_prct : 0.878040662678628"
## [1] "Variables auto-corrélées à plus de 0.7 : pfnb_prct et bfnb_prct : 0.734066280198133"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et dgf_prct : 0.946561903962206"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et Dotation_forfaitaireN_prct : 0.991748574148787"
## [1] "Variables auto-corrélées à plus de 0.7 : DotationNmoins1PerimN_prct et DSR_pereq_global_prct : 0.775196322883291"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et dgf_prct : 0.954457043664577"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DotationNmoins1PerimN_prct : 0.991748574148787"
## [1] "Variables auto-corrélées à plus de 0.7 : Dotation_forfaitaireN_prct et DSR_pereq_global_prct : 0.781293401890802"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_major_prct : 0.793043813747828"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_princip_prct et DNP_tot_prct : 0.979764635015164"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_princip_prct : 0.793043813747828"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_major_prct et DNP_tot_prct : 0.897271624224915"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_princip_prct : 0.979764635015164"
## [1] "Variables auto-corrélées à plus de 0.7 : DNP_tot_prct et DNP_major_prct : 0.897271624224915"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et dgf_prct : 0.795902189919594"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et DotationNmoins1PerimN_prct : 0.775196322883291"
## [1] "Variables auto-corrélées à plus de 0.7 : DSR_pereq_global_prct et Dotation_forfaitaireN_prct : 0.781293401890802"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_RETR1564_RT et P09_POP6074Y_RT : 0.740932400425599"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_ACT1564_Agr_RT et C09_EMPLT_AGRI_RT : 0.714166544320671"
## [1] "Variables auto-corrélées à plus de 0.7 : C09_EMPLT_AGRI_RT et C09_ACT1564_Agr_RT : 0.714166544320671"
## [1] "Variables auto-corrélées à plus de 0.7 : P09_POP6074Y_RT et P09_RETR1564_RT : 0.740932400425599"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_IMP_NET_RT : 0.879492075658382"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_Rev_Fisc_RT et P11_FoyFisc_Imp_RT : 0.734321656244132"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_IMP_NET_RT et P11_Rev_Fisc_RT : 0.879492075658382"
## [1] "Variables auto-corrélées à plus de 0.7 : P11_FoyFisc_Imp_RT et P11_Rev_Fisc_RT : 0.734321656244132"

Les variables envisagées ne sont pas trop auto-corrélées entre elles si on se réfère aux seuils préconisés dans la littérature.

Si jamais des variables sont trop autocorrélées, on peut les supprimer (c’est ce qui est fait dans la section suivante) puis relancer une étude de l’auto-corrélation.

variables_correlations <- variables_correlations[variables_correlations!="impo1_prct"]
variables_correlations <- variables_correlations[variables_correlations!="Pop_INSEE_prct"]
variables_correlations <- variables_correlations[variables_correlations!="Pop_DGF_prct"]
variables_correlations <- variables_correlations[variables_correlations!="Dotation_forfaitaireN_prct"]
variables_correlations <- variables_correlations[variables_correlations!="pfb_prct"]
variables_correlations <- variables_correlations[variables_correlations!="DNP_princip_prct"]
variables_correlations <- variables_correlations[variables_correlations!="DNP_major_prct"]
variables_correlations <- variables_correlations[variables_correlations!="DNP_tot_prct"]
variables_correlations <- variables_correlations[variables_correlations!="DSR_pereq_global_prct"]
variables_correlations <- variables_correlations[variables_correlations!="P11_POT_FIN_RT"]
variables_correlations <- variables_correlations[variables_correlations!="P11_Rev_Fisc_RT"]
variables_correlations <- variables_correlations[variables_correlations!="P11_IMP_NET_RT"]
variables_correlations <- variables_correlations[variables_correlations!="res1_prct"]
variables_correlations <- variables_correlations[variables_correlations!="depinv_prct"]
variables_correlations <- variables_correlations[variables_correlations!="tth_prct"]
variables_correlations <- variables_correlations[variables_correlations!="prod_prct"]
variables_correlations <- variables_correlations[variables_correlations!="res2_prct"]
variables_correlations <- variables_correlations[variables_correlations!="remb_prct"]
variables_correlations <- variables_correlations[variables_correlations!="pfnb_prct"]

rm(tmp, tmp2, res1, col, i, j, seuil)

4.2 Calcul de la Variance Inflaction Factor (VIF)

D’après la littérature [tenenhaus2007] , il faut calculer la Variance Inflation Factor (VIF) et supprimer les variables avec une VIF supérieure à 3 (voire supérieur à 2).

Pour la mise en œuvre, cf. cette page : https://larmarange.github.io/guide-R/analyses/multicolinearite.html

Ces éléments seront donc présentés lors de la mise en œuvre de la régression.

5 Comparaison des valeurs des communes fusionnantes vis-à-vis des valeurs des communes inchangées

# Calculs des valeurs moyennes (pour comparaison)
# On identifie les variables en stock
variables_stock <- stringr::str_replace(variables_RT, "Y_RT", "")
variables_stock <- stringr::str_replace(variables_stock, "_RT", "")
variables_stock <- variables_stock[-which(variables_stock == "C09_EMP_CONC")]
variables_stock <- variables_stock[1:(length(variables_stock)-5)]

# Import données totales
# load("data/refdata.Rdata")
# rm(df_new)
# Sélection données utiles
pourmoyennesCAH <- na.omit(subset(df2011, COM_NOUV == "OUI"))
moyennesCAH <- data.frame()
i <- variables_stock[2]
ratio <- as.data.frame(read_excel("data-raw/meta.xlsx", sheet = "ratios"))
row.names(ratio) <- ratio$Numerator_Code

for (i in variables_stock){
  a <- ratio[i, "Denominator_Code"]
  b <- ratio[i, "Coeff"]
  c <- sum(df2011[, i], na.rm = TRUE)
  d <- sum(df2011[, a], na.rm = TRUE)
  moyfr <- round(b * c/d, 2)
  e <- sum(pourmoyennesCAH[, i], na.rm = TRUE)
  f <- sum(pourmoyennesCAH[, a], na.rm = TRUE)
  moyCFus <- round(b * e/f, 2)
  g <- ratio[i, "CODE"]
  h <- c(g, moyfr, moyCFus)
  moyennesCAH <- rbind(moyennesCAH, h, stringsAsFactors= FALSE)
  rm( a, b, c, d, e, f, g, h, moyfr, moyCFus)
}

colnames(moyennesCAH) <- c("Variable", "France", "CommunesFusionnantes")
moyennesCAH$Variable <- as.factor(moyennesCAH$Variable)
moyennesCAH$France <- as.numeric(moyennesCAH$France)
moyennesCAH$CommunesFusionnantes <- as.numeric(moyennesCAH$CommunesFusionnantes)
moyennesCAH$DiffComFusComFr <- moyennesCAH$CommunesFusionnantes - moyennesCAH$France

aa<- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = France), stat = "identity") + coord_flip()
bb <- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = CommunesFusionnantes), stat = "identity") + coord_flip()
cc <- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = DiffComFusComFr), stat = "identity") + coord_flip() + ylab(paste0("Différence Communes fusionnantes - ", "France"))
cowplot::plot_grid(aa, bb, cc, ncol = 1, nrow = 3)

# Pour variables politiques
# Calculs des valeurs moyennes (pour comparaison)
# On identifie les variables en stock
VarCAHBrutes <- stringr::str_replace(variables_elect, "prct_insc", "nbre_voix")
# Quelques cas particuliers
VarCAHBrutes <- stringr::str_replace(VarCAHBrutes, "Abst_nbre_voix", "Abst_nbre")
VarCAHBrutes <- stringr::str_replace(VarCAHBrutes, "BlcsNuls_nbre_voix", "BlcsNuls_nbre")

elect <- substr(VarCAHBrutes[1], start = 1, stop = 9)

# Sélection données utiles
pourmoyennesCAH <- na.omit(subset(df2011, COM_NOUV == "OUI"))
somme_inscr_CFus <- sum(pourmoyennesCAH[, paste0(elect, "_Inscr_nbre")], na.rm = TRUE)
somme_inscr_ttes_com <- sum(na.omit(df2011[, paste0(elect, "_Inscr_nbre")]))

moyennesCAH <- data.frame()

i <- VarCAHBrutes[3]

for (i in VarCAHBrutes){
  somme_ttes_com <- sum(df2011[, i], na.rm = TRUE)
  moy_ttes_com <- round(100*somme_ttes_com/somme_inscr_ttes_com, 2)
  somme_CFus <- sum(pourmoyennesCAH[, i], na.rm = TRUE)
  moy_Cfus <- round(100*somme_CFus/somme_inscr_CFus, 2)
  variable <- variables_elect[ which(VarCAHBrutes == i) ]
  ligne <- c(variable, moy_ttes_com, moy_Cfus)
  moyennesCAH <- rbind(moyennesCAH, ligne, stringsAsFactors= FALSE)
  
}
rm (somme_inscr_CFus, somme_inscr_ttes_com, somme_ttes_com, moy_ttes_com, somme_CFus, moy_Cfus, ligne, variable)

colnames(moyennesCAH) <- c("Variable", "Ensemble_étudié", "CommunesFusionnantes")
moyennesCAH$Variable <- as.factor(moyennesCAH$Variable)
moyennesCAH$Ensemble_étudié <- as.numeric(moyennesCAH$Ensemble_étudié)
moyennesCAH$CommunesFusionnantes <- as.numeric(moyennesCAH$CommunesFusionnantes)
moyennesCAH$DiffComFusComFr <- moyennesCAH$CommunesFusionnantes - moyennesCAH$Ensemble_étudié

aa<- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = Ensemble_étudié), stat = "identity") + coord_flip() + ylab("France")
bb <- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = CommunesFusionnantes), stat = "identity") + coord_flip() + ylab("Communes fusionnantes")
cc <- ggplot(data = moyennesCAH) +
  geom_bar(aes(x = Variable, y = DiffComFusComFr), stat = "identity") + coord_flip() + ylab(paste0("Différence Communes fusionnantes - ", "France"))
cowplot::plot_grid(aa, bb, cc, ncol = 1, nrow = 3)

6 Régression logistique classique

Cf. par exemple https://larmarange.github.io/analyse-R/regression-logistique.html.

On fait le choix d’analyser ces données sous l’angle d’une régression logistique. Il s’agit d’observer si des variables augmentent ou diminuent la probabilité d’avoir une situation (fusion ou non). (Les tests de régression logistique sont utilisées parfois pour observer si des caractéristiques, qui seraient donc alors des comorbidités, augmentent ou non la probabilité d’avoir un décès.)

6.1 Variables explicatives envisagées

Les variables explicatives envisagées sont celles qui ont été pointées comme caractérisant davantage les communes nouvelles lors des études univariées, bivariées, ou lors de catégorisations [@bideau2019, @bideau2022c].

Le choix des variables explicatives envisagées est déterminant, comme cela a pu être montré par exemple par @afsa2016. En effet, le modèle logit présente la probabilité d’estimer la probabilité qu’un évènement se produise (dans notre cas, fusion ou non) mais il seulement en prenant en compte les données présentes dans le modèle. Dans le cas présenté par @afsa2016, si on n’intègre pas les CSP des parents, le fait d’être en ZEP a un impact très négatif sur l’orientation en seconde. Si on intègre ces CSP, la significativité de la ZEP disparaît presque intégralement, si on ajoute le niveau en 6e, le fait d’être en ZEP redevient significatif mais un facteur favorable pour cette orientation.

# Si l'on souhaite avoir des noms de variables plus explicites, il faut ajouter des étiquettes des variables avec var_label de l'extension labelled. Par exemple :
# var_label(d$sport) <- "Pratique du sport ?"

6.2 Modalités de la variable d’intérêt

Il est d’abord pertinent de vérifier que ce qui va être perçu, dans notre variable d’intérêt (ici la fusion ou non) comme la modalité de référence est bien la situation normale, la plus fréquente. Dans notre cas, la très grandes majorité des communes françaises n’a pas connu de fusion, c’est donc la non-fusion qui est notre modalité de référence.

# Nécessité de passer la variable d'intérêt en facteur pour la fonction glm suivante
df_reg$COM_NOUV <- as.factor(df_reg$COM_NOUV)
class(df_reg$COM_NOUV)

## [1] "factor"

freq(df_reg$COM_NOUV)

##         n    % val%
## NON 33537 92.9 92.9
## OUI  2571  7.1  7.1

# Vérification que la variable la plus fréquente est bien la première (doit être la variable de référence)
freq(df_reg$COM_NOUV, sort = "dec")

##         n    % val%
## NON 33537 92.9 92.9
## OUI  2571  7.1  7.1

row.names(freq(df_reg$COM_NOUV, sort = "dec"))

## [1] "NON" "OUI"

ordre <- row.names(freq(df_reg$COM_NOUV, sort = "dec"))
df_reg$COM_NOUV <- factor (df_reg$COM_NOUV, ordre)

# Si on veut forcer dans un sens :
# df_reg$COM_NOUV <- relevel (df_reg$COM_NOUV, "NON", "OUI")

6.3 Mise en œuvre de la régression logistique

Choix des variables qui doivent éviter celles pointées comme auto-corrélées plus haut. On affiche également le VIF pour limiter au maximum les biais d’interprétation.

# Les variables explicatives sont mentionnées juste après la variable d'intérêt (ici, "COM_NOUV")
reg <- glm(COM_NOUV ~ 
             # dgf_prct + perso_prct + equip_prct + dette_prct,
             # dgf_prct + perso_prct + equip_prct,
             # P09_RETR1564_RT + C09_ACT1564_Agr_RT + P09_POP3044Y_RT + C09_ACTOCC_OUT_RT + P09_CHOM1564_RT + dgf_prct + perso_prct + charge_prct,
             # P09_RETR1564_RT + C09_ACT1564_Agr_RT + P09_POP3044Y_RT + C09_ACTOCC_OUT_RT + P09_CHOM1564_RT + dgf_prct + perso_prct,
             # P09_RETR1564_RT + C09_ACT1564_Agr_RT + C09_ACTOCC_OUT_RT + P09_CHOM1564_RT,
             # P09_ETUD1564_RT + C09_ACT1564_Ouvr_RT + C09_ACTOCC_OUT_RT + P09_CHOM1564_RT,
             # C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT + dgf_prct,
             # C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT_quali + dgf_prct_quali, # Effets faibles mais significatifs
             # P09_ETUD1564_RT_quali + C09_ACT1564_Ouvr_RT_quali + C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT_quali,
             # C09_ACTOCC_OUT_RT_quali,
             # C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT_quali + P09_POP1529Y_RT_quali + P09_POP_quali + PR2012_T1_Abst_prct_insc_quali + PR2012_T1_prct_insc_LE.PEN_quali + CATAEU2010 + CODE_DEPT,
             P09_ETUD1564_RT_quali + C09_ACT1564_Ouvr_RT_quali + C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT_quali + P09_POP1529Y_RT_quali + P09_POP_quali + PR2012_T1_Abst_prct_insc_quali + PR2012_T1_prct_insc_LE.PEN_quali + dgf_prct_quali + perso_prct_quali + CATAEU2010,

           data = df_reg, family = binomial(logit), na.action = na.exclude)

summary(reg)

## 
## Call:
## glm(formula = COM_NOUV ~ P09_ETUD1564_RT_quali + C09_ACT1564_Ouvr_RT_quali + 
##     C09_ACTOCC_OUT_RT_quali + P09_CHOM1564_RT_quali + P09_POP1529Y_RT_quali + 
##     P09_POP_quali + PR2012_T1_Abst_prct_insc_quali + PR2012_T1_prct_insc_LE.PEN_quali + 
##     dgf_prct_quali + perso_prct_quali + CATAEU2010, family = binomial(logit), 
##     data = df_reg, na.action = na.exclude)
## 
## Coefficients:
##                                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                        -2.290702   0.128769 -17.789  < 2e-16 ***
## P09_ETUD1564_RT_qualiQ1             0.150951   0.069550   2.170 0.029976 *  
## P09_ETUD1564_RT_qualiQ2             0.149080   0.063980   2.330 0.019801 *  
## P09_ETUD1564_RT_qualiQ4            -0.017124   0.067759  -0.253 0.800483    
## P09_ETUD1564_RT_qualiQ5            -0.100474   0.073157  -1.373 0.169624    
## C09_ACT1564_Ouvr_RT_qualiQ1        -0.320170   0.076280  -4.197 2.70e-05 ***
## C09_ACT1564_Ouvr_RT_qualiQ2        -0.256018   0.069906  -3.662 0.000250 ***
## C09_ACT1564_Ouvr_RT_qualiQ4         0.173714   0.063328   2.743 0.006086 ** 
## C09_ACT1564_Ouvr_RT_qualiQ5         0.247586   0.065060   3.806 0.000142 ***
## C09_ACTOCC_OUT_RT_qualiQ1          -0.021315   0.073515  -0.290 0.771858    
## C09_ACTOCC_OUT_RT_qualiQ2           0.011743   0.065364   0.180 0.857429    
## C09_ACTOCC_OUT_RT_qualiQ4          -0.146694   0.068674  -2.136 0.032672 *  
## C09_ACTOCC_OUT_RT_qualiQ5          -0.046489   0.068632  -0.677 0.498172    
## P09_CHOM1564_RT_qualiQ1            -0.029202   0.066519  -0.439 0.660661    
## P09_CHOM1564_RT_qualiQ2             0.019444   0.063797   0.305 0.760529    
## P09_CHOM1564_RT_qualiQ4            -0.085348   0.065241  -1.308 0.190802    
## P09_CHOM1564_RT_qualiQ5            -0.376529   0.072666  -5.182 2.20e-07 ***
## P09_POP1529Y_RT_qualiQ1            -0.529325   0.073271  -7.224 5.04e-13 ***
## P09_POP1529Y_RT_qualiQ2            -0.348289   0.066055  -5.273 1.34e-07 ***
## P09_POP1529Y_RT_qualiQ4            -0.043979   0.064097  -0.686 0.492622    
##  [ getOption("max.print") est atteint -- 29 lignes omises ]
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 17789  on 34464  degrees of freedom
## Residual deviance: 17004  on 34416  degrees of freedom
##   (1643 observations effacées parce que manquantes)
## AIC: 17102
## 
## Number of Fisher Scoring iterations: 6

# exp(coef(reg))
# exp(confint(reg))
# exp(cbind(coef(reg), confint(reg)))
# odds.ratio(reg)

# Pour vérifier les VIF
car::vif(reg)

##                                      GVIF Df GVIF^(1/(2*Df))
## P09_ETUD1564_RT_quali            1.325792  4        1.035880
## C09_ACT1564_Ouvr_RT_quali        1.262476  4        1.029563
## C09_ACTOCC_OUT_RT_quali          1.586738  4        1.059408
## P09_CHOM1564_RT_quali            1.217186  4        1.024872
## P09_POP1529Y_RT_quali            1.431604  4        1.045870
## P09_POP_quali                    2.111837  4        1.097950
## PR2012_T1_Abst_prct_insc_quali   1.182556  4        1.021181
## PR2012_T1_prct_insc_LE.PEN_quali 1.220717  4        1.025243
## dgf_prct_quali                   1.316906  4        1.035010
## perso_prct_quali                 1.576806  4        1.058577
## CATAEU2010                       2.351000  8        1.054881

which(car::vif(reg) > 1.9)

##  [1]  6 11 12 13 14 15 16 17 18 19 20 21 22

# Présentation d'un tableau plus propre
# tbl_regression(reg, exponentiate = TRUE)
add_vif(tbl_regression(reg, exponentiate = TRUE)) # Avec affichage des VIF

Caractéristique	OR¹	95% IC¹	p-valeur	GVIF¹	Adjusted GVIF^2,1
P09_ETUD1564_RT_quali				1,3	1,0
Q3	—	—
Q1	1,16	1,01 – 1,33	0,030
Q2	1,16	1,02 – 1,32	0,020
Q4	0,98	0,86 – 1,12	0,8
Q5	0,90	0,78 – 1,04	0,2
C09_ACT1564_Ouvr_RT_quali				1,3	1,0
Q3	—	—
Q1	0,73	0,62 – 0,84	<0,001
Q2	0,77	0,67 – 0,89	<0,001
Q4	1,19	1,05 – 1,35	0,006
Q5	1,28	1,13 – 1,46	<0,001
C09_ACTOCC_OUT_RT_quali				1,6	1,1
Q3	—	—
Q1	0,98	0,85 – 1,13	0,8
Q2	1,01	0,89 – 1,15	0,9
Q4	0,86	0,75 – 0,99	0,033
Q5	0,95	0,83 – 1,09	0,5
P09_CHOM1564_RT_quali				1,2	1,0
Q3	—	—
Q1	0,97	0,85 – 1,11	0,7
Q2	1,02	0,90 – 1,16	0,8
Q4	0,92	0,81 – 1,04	0,2
Q5	0,69	0,59 – 0,79	<0,001
P09_POP1529Y_RT_quali				1,4	1,0
Q3	—	—
Q1	0,59	0,51 – 0,68	<0,001
Q2	0,71	0,62 – 0,80	<0,001
Q4	0,96	0,84 – 1,08	0,5
Q5	1,09	0,95 – 1,24	0,2
P09_POP_quali				2,1	1,1
Q3	—	—
Q1	0,74	0,64 – 0,86	<0,001
Q2	0,88	0,78 – 1,00	0,049
Q4	0,92	0,81 – 1,04	0,2
Q5	0,94	0,81 – 1,09	0,4
PR2012_T1_Abst_prct_insc_quali				1,2	1,0
Q3	—	—
Q1	1,16	1,02 – 1,33	0,024
Q2	1,05	0,92 – 1,20	0,5
Q4	1,02	0,89 – 1,16	0,8
Q5	1,00	0,87 – 1,15	>0,9
PR2012_T1_prct_insc_LE.PEN_quali				1,2	1,0
Q3	—	—
Q1	1,32	1,16 – 1,50	<0,001
Q2	1,24	1,09 – 1,40	<0,001
Q4	0,83	0,73 – 0,94	0,004
Q5	0,45	0,38 – 0,52	<0,001
dgf_prct_quali				1,3	1,0
Q3	—	—
Q1	0,81	0,70 – 0,94	0,006
Q2	0,87	0,76 – 0,99	0,042
Q4	1,10	0,97 – 1,25	0,2
Q5	1,49	1,31 – 1,70	<0,001
perso_prct_quali				1,6	1,1
Q3	—	—
Q1	1,09	0,95 – 1,26	0,2
Q2	1,06	0,93 – 1,20	0,4
Q4	0,96	0,84 – 1,09	0,5
Q5	0,75	0,65 – 0,86	<0,001
CATAEU2010				2,4	1,1
112	—	—
400	1,28	1,11 – 1,47	<0,001
300	1,20	1,06 – 1,36	0,005
120	1,18	1,02 – 1,36	0,027
111	0,45	0,34 – 0,58	<0,001
221	1,66	1,29 – 2,11	<0,001
212	1,29	0,99 – 1,66	0,053
222	0,79	0,53 – 1,15	0,2
211	0,84	0,53 – 1,26	0,4
¹ OR = rapport de cotes, IC = intervalle de confiance, GVIF = Generalized Variance Inflation Factor
² GVIF^[1/(2*df)]

# Représentations graphiques du modèle
# ggcoef_model(reg, exponentiate = FALSE)
forest_model(reg)

# Représentation graphique des effets
# plot(allEffects(reg))

plot_grid(plotlist = plot(ggeffect(reg)))

# Pour n'afficher qu'une variable :
# plot(ggeffect(reg, "_nom_de_la_variable"))

screenreg(reg, stars = c(0.01, 0.05, 0.1))

## 
## ================================================
##                                     Model 1     
## ------------------------------------------------
## (Intercept)                            -2.29 ***
##                                        (0.13)   
## P09_ETUD1564_RT_qualiQ1                 0.15 ** 
##                                        (0.07)   
## P09_ETUD1564_RT_qualiQ2                 0.15 ** 
##                                        (0.06)   
## P09_ETUD1564_RT_qualiQ4                -0.02    
##                                        (0.07)   
## P09_ETUD1564_RT_qualiQ5                -0.10    
##                                        (0.07)   
## C09_ACT1564_Ouvr_RT_qualiQ1            -0.32 ***
##                                        (0.08)   
## C09_ACT1564_Ouvr_RT_qualiQ2            -0.26 ***
##                                        (0.07)   
## C09_ACT1564_Ouvr_RT_qualiQ4             0.17 ***
##                                        (0.06)   
## C09_ACT1564_Ouvr_RT_qualiQ5             0.25 ***
##                                        (0.07)   
## C09_ACTOCC_OUT_RT_qualiQ1              -0.02    
##                                        (0.07)   
## C09_ACTOCC_OUT_RT_qualiQ2               0.01    
##                                        (0.07)   
## C09_ACTOCC_OUT_RT_qualiQ4              -0.15 ** 
##                                        (0.07)   
## C09_ACTOCC_OUT_RT_qualiQ5              -0.05    
##                                        (0.07)   
## P09_CHOM1564_RT_qualiQ1                -0.03    
##                                        (0.07)   
## P09_CHOM1564_RT_qualiQ2                 0.02    
##                                        (0.06)   
## P09_CHOM1564_RT_qualiQ4                -0.09    
##                                        (0.07)   
## P09_CHOM1564_RT_qualiQ5                -0.38 ***
##                                        (0.07)   
## P09_POP1529Y_RT_qualiQ1                -0.53 ***
##                                        (0.07)   
## P09_POP1529Y_RT_qualiQ2                -0.35 ***
##                                        (0.07)   
## P09_POP1529Y_RT_qualiQ4                -0.04    
##                                        (0.06)   
## P09_POP1529Y_RT_qualiQ5                 0.08    
##                                        (0.07)   
## P09_POP_qualiQ1                        -0.30 ***
##                                        (0.08)   
## P09_POP_qualiQ2                        -0.13 ** 
##                                        (0.06)   
## P09_POP_qualiQ4                        -0.09    
##                                        (0.06)   
## P09_POP_qualiQ5                        -0.06    
##                                        (0.08)   
## PR2012_T1_Abst_prct_insc_qualiQ1        0.15 ** 
##                                        (0.07)   
## PR2012_T1_Abst_prct_insc_qualiQ2        0.05    
##                                        (0.07)   
## PR2012_T1_Abst_prct_insc_qualiQ4        0.02    
##                                        (0.07)   
## PR2012_T1_Abst_prct_insc_qualiQ5        0.00    
##                                        (0.07)   
## PR2012_T1_prct_insc_LE.PEN_qualiQ1      0.28 ***
##                                        (0.07)   
## PR2012_T1_prct_insc_LE.PEN_qualiQ2      0.21 ***
##                                        (0.06)   
## PR2012_T1_prct_insc_LE.PEN_qualiQ4     -0.19 ***
##                                        (0.07)   
## PR2012_T1_prct_insc_LE.PEN_qualiQ5     -0.81 ***
##                                        (0.08)   
## dgf_prct_qualiQ1                       -0.21 ***
##                                        (0.07)   
## dgf_prct_qualiQ2                       -0.14 ** 
##                                        (0.07)   
## dgf_prct_qualiQ4                        0.09    
##                                        (0.07)   
## dgf_prct_qualiQ5                        0.40 ***
##                                        (0.07)   
## perso_prct_qualiQ1                      0.09    
##                                        (0.07)   
## perso_prct_qualiQ2                      0.06    
##                                        (0.07)   
## perso_prct_qualiQ4                     -0.04    
##                                        (0.06)   
## perso_prct_qualiQ5                     -0.29 ***
##                                        (0.07)   
## CATAEU2010400                           0.24 ***
##                                        (0.07)   
## CATAEU2010300                           0.18 ***
##                                        (0.06)   
## CATAEU2010120                           0.16 ** 
##                                        (0.07)   
## CATAEU2010111                          -0.80 ***
##                                        (0.13)   
## CATAEU2010221                           0.51 ***
##                                        (0.13)   
## CATAEU2010212                           0.26 *  
##                                        (0.13)   
## CATAEU2010222                          -0.23    
##                                        (0.20)   
## CATAEU2010211                          -0.18    
##                                        (0.22)   
## ------------------------------------------------
## AIC                                 17102.05    
## BIC                                 17515.99    
## Log Likelihood                      -8502.03    
## Deviance                            17004.05    
## Num. obs.                           34465       
## ================================================
## *** p < 0.01; ** p < 0.05; * p < 0.1

# texreg(reg, stars = c(0.01, 0.05, 0.1)) # Pour export LaTeX

6.4 [À vérifier] Lecture du modèle

Les odds ratio (littéralement rapport des cotes, ou rapport des chances ou rapport des risques relatiffs) permettent de décrire la significativité pratique. Celle-ci désigne l’impact d’une variable. Il ne faut pas la confondre avec la significativité statistique, qui décrit « le degré de certitude avec lequel on peut affirmer qu’une variable influe » (@afsa2016 p. 39).

À noter que, par exemple pour @afsa2016, l’expression « toutes choses égales par ailleurs » est à éviter. En effet, ce n’est le cas que si les variables étudiées permettent de décrire parfaitement le phénomène. « les résultats des estimations sont conditionnels à la liste des variables x introduites dans le modèle, c’est-à-dire qu’ils dépendent des variables introduites » (@afsa2016 p. 53). Si des variables, non mesurées, influent à la fois sur les variables du modèle et la variable étudiée, cela brouille l’analyse. Ainsi, pour les analyses présentées ici, c’est plutôt l’expression “Toutes choses égales quant aux autres variables introduites” qui est à utiliser.

Un Odds ratio égal à 1 signifie que la variable d’intérêt (la fusion ici) est indépendante de la variable explicative envisagée. Si l’OR est supérieur à 1, une commune ayant une variable explicative plus élevée ont plus de chance de fusionner, s’il est inférieur à 1, les communes avec une variable explicative plus faible ont moins de chance de fusionner.

6.5 Sauvegarde données modèle

try(save(df2011, df_new, df_reg, reg, file = "data/refdata_modele.Rdata"))

load("data/refdata_modele.Rdata")

7 [À faire] Étude des résidus (y compris avec indice de Moran pour une auto-corrélation spatiale des résidus ?)

Non-stationnarité spatiale : lorsque les paramètres permettant d’expliquer un phénomènes sont différents en fonction des espaces. Bonne figure p. 185 (Feuillet et al, 2019) où en regardant deux groupes distincts on aperçoit un type de relation alors que lorsqu’on regarde les deux en même temps, on imagine une autre relation statistique. C’est le paradoxe de Simpson (qui n’est pas forcément situé).

Il est donc important voire intéressant de cartographier ou de représenter graphiquement les résidus. Cela permet potentiellement de valider ou invalider le modèle mais également de réfléchir aux espaces où ce modèle fonctionne le mieux ou le moins bien.

7.1 Étude résidus régression linéaire (peu pertinent dans le cas d’un logit)

Cf. R et Espace section 5.2.3 et particulièrement p. 91 sq pour aperçu des résidus : « La fonction plot() appliquée aux résultats d’un modèle de régression linéaire obtenus par la fonction lm() permet de représenter les quatre principales hypothèses au cœur de ce modèle :

— la normalité des résidus par rapport aux valeurs prédites (en haut à gauche de l’image)

— la normalité globale des résidus (en haut à droite de l’image)

— la corrélation entre les valeurs de la variable explicative et le carré des résidus standardisés (en bas à gauche de l’image)

— l’existence de valeurs extrêmes altérant l’estimation des paramètres (en bas à droite de l’image) »

Autre source très utilisée ci-dessous : https://sites.google.com/site/rgraphiques/4--stat/machine-learning-biostatistiques-analyse-de-donn%C3%A9es/r%C3%A9gressions-lin%C3%A9aires-avec-r?pli=1

# On veut éviter les phénomènes d'auto-corrélation : les résidus augmenteraient ensemble dans une zone donnée et seraient liés les uns aux résidus précédents ou suivants (cas typiques lorsque X est un enregistrement du temps). Vérification par le test de Durbin-Watson (indépendance : p-value > 0,05)
durbinWatsonTest(reg)
dwtest(reg)

# Visualisation des résidus
par(mfrow = c(2,2))
plot(reg)
par(mfrow = c(1,1))

# Les résidus doivent être répartis approximativement selon une courbe de Gauss
hist(residuals(reg),col="yellow",freq=F)

# Il faut, pour que la régression soit pertinente, observer un alignement entre les quantiles des résidus et les quantiles théoriques, donc globalement une ligne droite
plot(reg, which = 2)
# On peut vérifier avec le test de normalité de Shapiro-Wilk (normalité si p-value supérieure à 0,05)
shapiro.test(residuals(reg))


# Vérification que la variance des résidus est constante (distribution homogène)
# Pour que le modèle soit pertinent, il faut une courbe rouge plane et que les valeurs ne se regroupent pas.
plot(reg, which = 3)
# On peut vérifier avec le test de Breush-Pagan (homogénéité si p-value > 0,05)
# ncvTest(reg) # seulement pour objet lm
bptest(reg)
# Ou test de Goldfeld-Quandt (homogénéité si p-value > 0,05)
gqtest(reg)

plot(reg, which = 5)

7.2 Étude résidus logit

Cf. https://pmarchand1.github.io/ECL7102/notes_cours/9-Regression_logistique.pdf#page=9

binnedplot(fitted(reg), residuals(reg, type = "response"))

7.3 Tentative cartographie résidus

# residuals (reg)
# On intègre les résidus aux données
df_reg$residus <- residuals (reg)
par(mar=c(0,0,1.2,0), mfrow = c(1,2))

plot(st_geometry(dep), col = NA)

choroLayer(x = df_reg , var = "residus", 
           method = "quantile", nclass = 6,
           col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
           border = NA,
           colNA = "white",
           legend.pos = "topleft", legend.values.rnd = 2,
           legend.title.txt = "Résidus",
           add = TRUE)

layoutLayer(title = "Cartographie des résidus",
            author = "G. Bideau, 2024.",
            sources = "Sources : ")

typoLayer(x = df_reg, var = "COM_NOUV",
          legend.values.order = c("OUI", "NON"),
          legend.pos = "topleft",
          col=c("red","blue"),
          legend.title.txt = "Communes nouvelles",
          border = NA)

layoutLayer(title = "Communes nouvelles 2012-2024(01)",
            author = "G. Bideau, 2024.",
            sources = "Sources : ")

par(mar = c(0,0,0,0), mfrow = c(1,1))

7.4 Test de Moran sur auto-corrélation spatiale des résidus (package rgeoda)

test_geom <- df_reg

# test_geom$COM_NOUV2<- if_else(test_geom$COM_NOUV == "OUI", 1, 0)
# variable <- "COM_NOUV2"
variable <- "residus"
nom_variable <- "Résidus"

df_variable <- test_geom[variable]

# Définition de la matrice de contiguïté
queen_w <- queen_weights(test_geom, order=1, include_lower_order = FALSE, precision_threshold = 0)

lisa <- local_moran(queen_w, df_variable)

# To get the values of the local Moran:
lms <- lisa_values(gda_lisa = lisa)

# To get the pseudo-p values of significance of local Moran computation:
pvals <- lisa_pvalues(lisa)

# To get the cluster indicators of local Moran computation:
cats <- lisa_clusters(lisa, cutoff = 0.05)

# Labels
lbls <- lisa_labels(lisa)


test_geom$LISA <- lms
test_geom$pvals <- pvals
test_geom$cats <- cats

# On rend les catégories plus lisibles en remplaçant le numéro par l'intitulé
test_geom$cats <- as.factor(test_geom$cats)
num_levels <- as.numeric(levels(test_geom$cats))
levels(test_geom$cats) <- lbls[num_levels + 1]
print(levels(test_geom$cats))

## [1] "High-High" "Isolated"

test_geom$LISA_retenu <- ifelse(test_geom$pvals <= 0.1, yes = test_geom$LISA, no = NA)

par(mar = c(0,0,0,0), mfrow = c(1,2))

couleurs <- c("#f7f7f7", "#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#ffff33")
# Cf. https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6 pour composer les palettes
typoLayer(x = test_geom , var = "cats",
           # col = carto.pal(pal1 = "blue.pal", pal2 = "red.pal", n1 = 3, n2 = 3),
           col = couleurs[1:length(levels(test_geom$cats))],
           border = NA,
           legend.pos = "topleft",
           legend.values.order = levels(test_geom$cats),
           legend.title.txt = paste0("Indice de Moran locaux significatifs\nsur la variable « ", nom_variable, " »\n(pvalue < 0,1)"))


typoLayer(x = test_geom, var = "COM_NOUV",
          legend.values.order = c("OUI", "NON"),
          col=c("red","blue"),
          border = NA)

# plot(test_geomCN$geometry, col = NA, border = "black", lwd = 1, add = TRUE)

par(mar = c(0,0,0,0), mfrow = c(1,1))

rm(variable, df_variable, lms, pvals, cats, lbls, num_levels, couleurs)

8 [À vérifier, tester] Forward stepwise ou backward stepwise

Possibilité de sélectionner automatiquement les variables à choisir pour une régression logistique [feuillet2019] :

régression pas à pas ascendante (on ajoute des variables à partir d’un tableau de données) : forward stepwise ;
régression pas à pas descendante (on enlève) : backward stepwise.

Guides pour mettre en place ces méthodes : https://www.statology.org/stepwise-regression-r/ pour une présentation un peu rapide ou (déroulé suivi ici) : https://rstudio-pubs-static.s3.amazonaws.com/205694_3b195f29e9504d23aeb483ff1ffafeba.html, en particulier à partir du point 3.4.

8.1 Sélection des données

# Sélection d'une partie de df_reg pour le step by step

variables_sbs <- c(colnames(df_reg[which(str_detect(colnames(df_reg), "_prct") == TRUE)]),
                            colnames(df_reg[which(str_detect(colnames(df_reg), "_RT") == TRUE)]))
variables_sbs <- c(colnames(df_reg[which(str_detect(colnames(df_reg), "_prct") == TRUE)]),
                            colnames(df_reg[which(str_detect(colnames(df_reg), "_RT") == TRUE)]), "CATAEU2010", "CODE_DEPT")
variables_sbs <- c(variables_sbs, variables_elect)
variables_sbs <- variables_sbs[variables_sbs!= "geometry"] # On supprime la variable "geometry" qui est intégrée sans que cela soit voulu

# On ne conserve que les variables qu'on a validé du point de vue du risque d'auto-corrélation
variables_sbs <- intersect(variables_sbs, variables_correlations) 

df_reg_sbs <- df_reg[, c("COM_NOUV", variables_sbs)]

df_reg_sbs <- na.omit(df_reg_sbs)

# Suppression des géométries
df_reg_sbs <- st_set_geometry (df_reg_sbs, NULL)

8.2 Régression pas à pas

# Pour procéder avec la sélection ascendante (Forward selection) on doit d’abord spécifier le modèle de départ
m0 <- glm(COM_NOUV ~ 1, data=df_reg_sbs, family = binomial(logit), na.action = na.exclude)  # choix du modèle avec constante seulement
# m0

# Pour utiliser la commande step on doit spécifier le modèle de départ m0 et le modèle maximal mf qui peut être le modèle complet.
mf <- glm(COM_NOUV ~ ., data=df_reg_sbs, family = binomial(logit), na.action = na.exclude)  # choix du modèle avec constante seulement
# mf

# La sélection ascendante utilisant le critère AIC
# step(m0, scope=list(lower=m0, upper=mf), data = df_reg_sbs, direction="forward")
reg_sbs_forward <- step(m0, scope=list(lower=m0, upper=mf), data = df_reg_sbs, direction="forward", trace = FALSE)
add_vif(tbl_regression(reg_sbs_forward, exponentiate = TRUE)) # Pour afficher les VIF

Caractéristique	OR¹	95% IC¹	p-valeur	VIF¹
PR2012_T1_prct_insc_BAYROU	1,01	0,99 – 1,03	0,5	3,0
C09_ACT1564_Ouvr_RT	1,01	1,01 – 1,02	<0,001	1,5
PR2012_T1_prct_insc_LE.PEN	0,92	0,90 – 0,93	<0,001	4,9
bfnb_prct	1,01	1,00 – 1,01	<0,001	1,9
P09_POP4559Y_RT	0,99	0,98 – 1,01	0,5	1,8
PR2012_T1_prct_insc_SARKOZY	0,99	0,97 – 1,00	0,15	7,5
dgf_prct	1,02	1,02 – 1,03	<0,001	2,3
achat_prct	0,97	0,97 – 0,98	<0,001	1,6
fdr_prct	1,00	1,00 – 1,00	<0,001	1,7
PotFin4taxes_prct	1,00	1,00 – 1,01	<0,001	2,2
pth_prct	0,97	0,96 – 0,98	<0,001	1,8
P09_POP0014Y_RT	1,06	1,04 – 1,08	<0,001	3,6
DSU_tot_prct	0,78	0,67 – 0,89	<0,001	1,1
subv_prct	1,02	1,01 – 1,03	0,002	1,1
P09_POP1529Y_RT	1,06	1,04 – 1,08	<0,001	2,2
P09_ETUD1564_RT	0,98	0,97 – 0,99	<0,001	1,3
P09_RETR1564_RT	1,02	1,01 – 1,03	<0,001	1,8
equip_prct	0,99	0,99 – 1,0	<0,001	2,5
perso_prct	0,99	0,98 – 0,99	<0,001	1,7
C09_ACT1564_Agr_RT	0,99	0,98 – 1,00	<0,001	1,6
P09_POP75PY_RT	1,03	1,01 – 1,04	0,001	3,0
C09_EMPLT_CTS_RT	1,00	1,00 – 1,00	0,003	1,2
C09_ACT1564_Cadr_RT	0,99	0,98 – 1,00	0,007	1,4
C09_ACT1564_ProfInt_RT	0,99	0,99 – 1,00	0,006	1,4
fctva_prct	0,99	0,98 – 1,00	0,039	1,1
P09_CHOM1564_RT	0,99	0,97 – 1,00	0,025	1,3
impo2_prct	0,99	0,99 – 1,00	0,036	1,2
PR2012_T1_prct_insc_HOLLANDE	0,96	0,95 – 0,98	<0,001	6,6
PR2012_T1_prct_insc_MÉLENCHON	0,96	0,94 – 0,98	<0,001	3,5
PR2012_T1_Abst_prct_insc	0,96	0,95 – 0,98	<0,001	3,3
annu_prct	0,99	0,98 – 1,00	0,010	2,0
dette_prct	1,00	1,00 – 1,00	0,072	1,7
P11_DGF_RT	1,00	1,00 – 1,00	0,019	1,4
bth_prct	1,00	1,00 – 1,00	0,070	3,0
C09_EMPLT_INDUS_RT	1,00	1,00 – 1,00	0,082	1,2
¹ OR = rapport de cotes, IC = intervalle de confiance, VIF = Variance Inflation Factor

# La sélection descendante
# step(mf, data = df_reg_sbs, direction="backward")
reg_sbs_backward <- step(mf, data = df_reg_sbs, direction="backward", na.action = na.omit, trace = FALSE)
add_vif(tbl_regression(reg_sbs_backward, exponentiate = TRUE)) # Pour afficher les VIF

Caractéristique	OR¹	95% IC¹	p-valeur	VIF¹
impo2_prct	0,99	0,99 – 1,00	0,036	1,2
dgf_prct	1,02	1,02 – 1,03	<0,001	2,3
perso_prct	0,99	0,98 – 0,99	<0,001	1,7
achat_prct	0,97	0,97 – 0,98	<0,001	1,6
subv_prct	1,02	1,01 – 1,03	0,001	1,1
fctva_prct	0,99	0,98 – 1,00	0,039	1,1
equip_prct	0,99	0,99 – 1,0	<0,001	2,5
dette_prct	1,00	1,00 – 1,00	0,071	1,7
annu_prct	0,99	0,98 – 1,00	0,011	2,0
fdr_prct	1,00	1,00 – 1,00	<0,001	1,7
bth_prct	1,00	1,00 – 1,00	0,070	3,1
bfnb_prct	1,01	1,00 – 1,01	<0,001	1,9
pth_prct	0,97	0,96 – 0,98	<0,001	1,8
PotFin4taxes_prct	1,00	1,00 – 1,01	<0,001	2,2
DSU_tot_prct	0,79	0,67 – 0,89	<0,001	1,1
PR2012_T1_Abst_prct_insc	0,96	0,95 – 0,97	<0,001	1,8
PR2012_T1_prct_insc_LE.PEN	0,91	0,90 – 0,92	<0,001	2,2
PR2012_T1_prct_insc_SARKOZY	0,98	0,97 – 1,00	0,008	3,9
PR2012_T1_prct_insc_MÉLENCHON	0,95	0,94 – 0,97	<0,001	2,1
PR2012_T1_prct_insc_HOLLANDE	0,96	0,95 – 0,97	<0,001	3,3
P09_CHOM1564_RT	0,99	0,97 – 1,00	0,024	1,3
P09_ETUD1564_RT	0,98	0,97 – 0,99	<0,001	1,2
P09_RETR1564_RT	1,02	1,01 – 1,03	<0,001	2,4
C09_ACT1564_Agr_RT	0,99	0,98 – 1,00	<0,001	1,6
C09_ACT1564_Cadr_RT	0,99	0,98 – 1,00	0,007	1,4
C09_ACT1564_ProfInt_RT	0,99	0,99 – 1,00	0,006	1,4
C09_ACT1564_Ouvr_RT	1,01	1,01 – 1,02	<0,001	1,5
C09_EMPLT_INDUS_RT	1,00	1,00 – 1,00	0,081	1,2
C09_EMPLT_CTS_RT	1,00	1,00 – 1,00	0,003	1,2
P09_POP3044Y_RT	0,94	0,93 – 0,96	<0,001	3,5
P09_POP4559Y_RT	0,94	0,92 – 0,95	<0,001	1,7
P09_POP6074Y_RT	0,94	0,93 – 0,96	<0,001	3,2
P09_POP75PY_RT	0,97	0,96 – 0,98	<0,001	2,2
P11_DGF_RT	1,00	1,00 – 1,00	0,020	1,4
¹ OR = rapport de cotes, IC = intervalle de confiance, VIF = Variance Inflation Factor

# La sélection pas a pas dans les deux sens
# step(m0, scope = list(upper=mf), data = df_reg_sbs,direction="both")
reg_sbs_both <- step(m0, scope = list(upper=mf), data = df_reg_sbs,direction="both", trace = FALSE)
add_vif(tbl_regression(reg_sbs_both, exponentiate = TRUE)) # Pour afficher les VIF

Caractéristique	OR¹	95% IC¹	p-valeur	VIF¹
C09_ACT1564_Ouvr_RT	1,01	1,01 – 1,02	<0,001	1,5
PR2012_T1_prct_insc_LE.PEN	0,91	0,90 – 0,92	<0,001	2,2
bfnb_prct	1,01	1,00 – 1,01	<0,001	1,9
PR2012_T1_prct_insc_SARKOZY	0,98	0,97 – 1,00	0,008	3,8
dgf_prct	1,02	1,02 – 1,03	<0,001	2,3
achat_prct	0,97	0,97 – 0,98	<0,001	1,6
fdr_prct	1,00	1,00 – 1,00	<0,001	1,7
PotFin4taxes_prct	1,00	1,00 – 1,01	<0,001	2,2
pth_prct	0,97	0,96 – 0,98	<0,001	1,8
P09_POP0014Y_RT	1,07	1,05 – 1,08	<0,001	2,2
DSU_tot_prct	0,79	0,67 – 0,89	<0,001	1,1
subv_prct	1,02	1,01 – 1,03	0,001	1,1
P09_POP1529Y_RT	1,06	1,05 – 1,08	<0,001	2,0
P09_ETUD1564_RT	0,98	0,97 – 0,99	<0,001	1,3
P09_RETR1564_RT	1,02	1,01 – 1,03	<0,001	1,7
equip_prct	0,99	0,99 – 1,0	<0,001	2,5
perso_prct	0,99	0,98 – 0,99	<0,001	1,7
C09_ACT1564_Agr_RT	0,99	0,98 – 1,00	<0,001	1,6
P09_POP75PY_RT	1,03	1,02 – 1,04	<0,001	2,3
C09_EMPLT_CTS_RT	1,00	1,00 – 1,00	0,003	1,2
C09_ACT1564_Cadr_RT	0,99	0,98 – 1,00	0,006	1,4
C09_ACT1564_ProfInt_RT	0,99	0,99 – 1,00	0,006	1,4
fctva_prct	0,99	0,98 – 1,00	0,040	1,1
P09_CHOM1564_RT	0,99	0,97 – 1,00	0,020	1,3
impo2_prct	0,99	0,99 – 1,00	0,036	1,2
PR2012_T1_prct_insc_HOLLANDE	0,96	0,95 – 0,97	<0,001	3,3
PR2012_T1_prct_insc_MÉLENCHON	0,95	0,94 – 0,97	<0,001	2,1
PR2012_T1_Abst_prct_insc	0,96	0,95 – 0,97	<0,001	1,8
annu_prct	0,99	0,98 – 1,00	0,010	2,0
dette_prct	1,00	1,00 – 1,00	0,069	1,7
P11_DGF_RT	1,00	1,00 – 1,00	0,020	1,4
bth_prct	1,00	1,00 – 1,00	0,070	3,0
C09_EMPLT_INDUS_RT	1,00	1,00 – 1,00	0,083	1,2
¹ OR = rapport de cotes, IC = intervalle de confiance, VIF = Variance Inflation Factor

# La sélection ascendante utilisant le F−test des modèles emboîtés.
# step(m0, scope=list(lower=m0, upper=mf), data = df_reg_sbs, direction="forward",test="F")
reg_sbs_forward_F <- step(m0, scope=list(lower=m0, upper=mf), data = df_reg_sbs, direction="forward",test="F", trace = FALSE)
add_vif(tbl_regression(reg_sbs_forward_F, exponentiate = TRUE)) # Pour afficher les VIF

Caractéristique	OR¹	95% IC¹	p-valeur	VIF¹
PR2012_T1_prct_insc_BAYROU	1,01	0,99 – 1,03	0,5	3,0
C09_ACT1564_Ouvr_RT	1,01	1,01 – 1,02	<0,001	1,5
PR2012_T1_prct_insc_LE.PEN	0,92	0,90 – 0,93	<0,001	4,9
bfnb_prct	1,01	1,00 – 1,01	<0,001	1,9
P09_POP4559Y_RT	0,99	0,98 – 1,01	0,5	1,8
PR2012_T1_prct_insc_SARKOZY	0,99	0,97 – 1,00	0,15	7,5
dgf_prct	1,02	1,02 – 1,03	<0,001	2,3
achat_prct	0,97	0,97 – 0,98	<0,001	1,6
fdr_prct	1,00	1,00 – 1,00	<0,001	1,7
PotFin4taxes_prct	1,00	1,00 – 1,01	<0,001	2,2
pth_prct	0,97	0,96 – 0,98	<0,001	1,8
P09_POP0014Y_RT	1,06	1,04 – 1,08	<0,001	3,6
DSU_tot_prct	0,78	0,67 – 0,89	<0,001	1,1
subv_prct	1,02	1,01 – 1,03	0,002	1,1
P09_POP1529Y_RT	1,06	1,04 – 1,08	<0,001	2,2
P09_ETUD1564_RT	0,98	0,97 – 0,99	<0,001	1,3
P09_RETR1564_RT	1,02	1,01 – 1,03	<0,001	1,8
equip_prct	0,99	0,99 – 1,0	<0,001	2,5
perso_prct	0,99	0,98 – 0,99	<0,001	1,7
C09_ACT1564_Agr_RT	0,99	0,98 – 1,00	<0,001	1,6
P09_POP75PY_RT	1,03	1,01 – 1,04	0,001	3,0
C09_EMPLT_CTS_RT	1,00	1,00 – 1,00	0,003	1,2
C09_ACT1564_Cadr_RT	0,99	0,98 – 1,00	0,007	1,4
C09_ACT1564_ProfInt_RT	0,99	0,99 – 1,00	0,006	1,4
fctva_prct	0,99	0,98 – 1,00	0,039	1,1
P09_CHOM1564_RT	0,99	0,97 – 1,00	0,025	1,3
impo2_prct	0,99	0,99 – 1,00	0,036	1,2
PR2012_T1_prct_insc_HOLLANDE	0,96	0,95 – 0,98	<0,001	6,6
PR2012_T1_prct_insc_MÉLENCHON	0,96	0,94 – 0,98	<0,001	3,5
PR2012_T1_Abst_prct_insc	0,96	0,95 – 0,98	<0,001	3,3
annu_prct	0,99	0,98 – 1,00	0,010	2,0
dette_prct	1,00	1,00 – 1,00	0,072	1,7
P11_DGF_RT	1,00	1,00 – 1,00	0,019	1,4
bth_prct	1,00	1,00 – 1,00	0,070	3,0
C09_EMPLT_INDUS_RT	1,00	1,00 – 1,00	0,082	1,2
¹ OR = rapport de cotes, IC = intervalle de confiance, VIF = Variance Inflation Factor

# La sélection descendante avec le F−test
# step(mf, data = df_reg_sbs, direction="backward")
reg_sbs_backward_F <- step(mf, data = df_reg_sbs, direction="backward", trace = FALSE)
add_vif(tbl_regression(reg_sbs_backward_F, exponentiate = TRUE)) # Pour afficher les VIF

Caractéristique	OR¹	95% IC¹	p-valeur	VIF¹
impo2_prct	0,99	0,99 – 1,00	0,036	1,2
dgf_prct	1,02	1,02 – 1,03	<0,001	2,3
perso_prct	0,99	0,98 – 0,99	<0,001	1,7
achat_prct	0,97	0,97 – 0,98	<0,001	1,6
subv_prct	1,02	1,01 – 1,03	0,001	1,1
fctva_prct	0,99	0,98 – 1,00	0,039	1,1
equip_prct	0,99	0,99 – 1,0	<0,001	2,5
dette_prct	1,00	1,00 – 1,00	0,071	1,7
annu_prct	0,99	0,98 – 1,00	0,011	2,0
fdr_prct	1,00	1,00 – 1,00	<0,001	1,7
bth_prct	1,00	1,00 – 1,00	0,070	3,1
bfnb_prct	1,01	1,00 – 1,01	<0,001	1,9
pth_prct	0,97	0,96 – 0,98	<0,001	1,8
PotFin4taxes_prct	1,00	1,00 – 1,01	<0,001	2,2
DSU_tot_prct	0,79	0,67 – 0,89	<0,001	1,1
PR2012_T1_Abst_prct_insc	0,96	0,95 – 0,97	<0,001	1,8
PR2012_T1_prct_insc_LE.PEN	0,91	0,90 – 0,92	<0,001	2,2
PR2012_T1_prct_insc_SARKOZY	0,98	0,97 – 1,00	0,008	3,9
PR2012_T1_prct_insc_MÉLENCHON	0,95	0,94 – 0,97	<0,001	2,1
PR2012_T1_prct_insc_HOLLANDE	0,96	0,95 – 0,97	<0,001	3,3
P09_CHOM1564_RT	0,99	0,97 – 1,00	0,024	1,3
P09_ETUD1564_RT	0,98	0,97 – 0,99	<0,001	1,2
P09_RETR1564_RT	1,02	1,01 – 1,03	<0,001	2,4
C09_ACT1564_Agr_RT	0,99	0,98 – 1,00	<0,001	1,6
C09_ACT1564_Cadr_RT	0,99	0,98 – 1,00	0,007	1,4
C09_ACT1564_ProfInt_RT	0,99	0,99 – 1,00	0,006	1,4
C09_ACT1564_Ouvr_RT	1,01	1,01 – 1,02	<0,001	1,5
C09_EMPLT_INDUS_RT	1,00	1,00 – 1,00	0,081	1,2
C09_EMPLT_CTS_RT	1,00	1,00 – 1,00	0,003	1,2
P09_POP3044Y_RT	0,94	0,93 – 0,96	<0,001	3,5
P09_POP4559Y_RT	0,94	0,92 – 0,95	<0,001	1,7
P09_POP6074Y_RT	0,94	0,93 – 0,96	<0,001	3,2
P09_POP75PY_RT	0,97	0,96 – 0,98	<0,001	2,2
P11_DGF_RT	1,00	1,00 – 1,00	0,020	1,4
¹ OR = rapport de cotes, IC = intervalle de confiance, VIF = Variance Inflation Factor

# La sélection pas à pas dans les deux sens avec le F−test
# step(m0, scope = list(upper=mf), data = df_reg_sbs, direction="both")
reg_sbs_both_F <- step(m0, scope = list(upper=mf), data = df_reg_sbs, direction="both", trace = FALSE)
add_vif(tbl_regression(reg_sbs_both_F, exponentiate = TRUE)) # Pour afficher les VIF

Caractéristique	OR¹	95% IC¹	p-valeur	VIF¹
C09_ACT1564_Ouvr_RT	1,01	1,01 – 1,02	<0,001	1,5
PR2012_T1_prct_insc_LE.PEN	0,91	0,90 – 0,92	<0,001	2,2
bfnb_prct	1,01	1,00 – 1,01	<0,001	1,9
PR2012_T1_prct_insc_SARKOZY	0,98	0,97 – 1,00	0,008	3,8
dgf_prct	1,02	1,02 – 1,03	<0,001	2,3
achat_prct	0,97	0,97 – 0,98	<0,001	1,6
fdr_prct	1,00	1,00 – 1,00	<0,001	1,7
PotFin4taxes_prct	1,00	1,00 – 1,01	<0,001	2,2
pth_prct	0,97	0,96 – 0,98	<0,001	1,8
P09_POP0014Y_RT	1,07	1,05 – 1,08	<0,001	2,2
DSU_tot_prct	0,79	0,67 – 0,89	<0,001	1,1
subv_prct	1,02	1,01 – 1,03	0,001	1,1
P09_POP1529Y_RT	1,06	1,05 – 1,08	<0,001	2,0
P09_ETUD1564_RT	0,98	0,97 – 0,99	<0,001	1,3
P09_RETR1564_RT	1,02	1,01 – 1,03	<0,001	1,7
equip_prct	0,99	0,99 – 1,0	<0,001	2,5
perso_prct	0,99	0,98 – 0,99	<0,001	1,7
C09_ACT1564_Agr_RT	0,99	0,98 – 1,00	<0,001	1,6
P09_POP75PY_RT	1,03	1,02 – 1,04	<0,001	2,3
C09_EMPLT_CTS_RT	1,00	1,00 – 1,00	0,003	1,2
C09_ACT1564_Cadr_RT	0,99	0,98 – 1,00	0,006	1,4
C09_ACT1564_ProfInt_RT	0,99	0,99 – 1,00	0,006	1,4
fctva_prct	0,99	0,98 – 1,00	0,040	1,1
P09_CHOM1564_RT	0,99	0,97 – 1,00	0,020	1,3
impo2_prct	0,99	0,99 – 1,00	0,036	1,2
PR2012_T1_prct_insc_HOLLANDE	0,96	0,95 – 0,97	<0,001	3,3
PR2012_T1_prct_insc_MÉLENCHON	0,95	0,94 – 0,97	<0,001	2,1
PR2012_T1_Abst_prct_insc	0,96	0,95 – 0,97	<0,001	1,8
annu_prct	0,99	0,98 – 1,00	0,010	2,0
dette_prct	1,00	1,00 – 1,00	0,069	1,7
P11_DGF_RT	1,00	1,00 – 1,00	0,020	1,4
bth_prct	1,00	1,00 – 1,00	0,070	3,0
C09_EMPLT_INDUS_RT	1,00	1,00 – 1,00	0,083	1,2
¹ OR = rapport de cotes, IC = intervalle de confiance, VIF = Variance Inflation Factor

8.3 Algorithme génétique

Section 3.5 de la page déjà citée : https://rstudio-pubs-static.s3.amazonaws.com/205694_3b195f29e9504d23aeb483ff1ffafeba.html#lalgorithme-genetique

Mais ne fonctionne pas : problème pour installer glmulti sur RStudio Server + message d’erreur “!Oversized candidate set.”

df_reg_gen <- df_reg_sbs[1:50,1:10]
res3 <- glmulti(COM_NOUV~.,data = df_reg_gen,level = 2,method = "g",fitfunction = lm,crit = 'aic',plotty = F)

#CN_regressions Les communes nouvelles françaises (2012-2024) : proposition d’un modèle explicatif à l’aide de régressions logistiques

Gabriel Bideau

19 juin 2024 23:04:51